maxvogt-analysis: Erlenbach

The code in this notebook is based on Martin Fleischmann’s 2021 workshop Capturing the Structure of Cities with Data Science (SDSC), licensed under the CC BY-SA 4.0 license. © 2021 Martin Fleischmann.

The geodata in this notebook is retrieved from OpenStreetMap. © OpenStreetMap Contributors.

This notebook contains geodata and analyses for the following items of the Max Vogt collection by Moritz Twente and Luisa Omonsky: - ERLB027

Open in an interactive in-browser environment:

Binder

Binder
import warnings

import geopandas
import libpysal
import momepy
import osmnx
import pandas

from clustergram import Clustergram

import matplotlib.pyplot as plt
from bokeh.io import output_notebook
from bokeh.plotting import show

Pick a place, ideally a town with a good coverage in OpenStreetMap and its local CRS.

place = 'Erlenbach ZH'
point = 47.30576, 8.59146
address = 'Bahnhofstrasse 28, 8703 Erlenbach ZH'
dist = 2000
local_crs = 'EPSG:2056'
geopandas.tools.geocode(address).explore()
Make this Notebook Trusted to load map: File -> Trust Notebook

Input data

Download data from OpenStreetMap.

Buildings

#buildings = osmnx.features.features_from_place(place, tags={'building':True})
#buildings = osmnx.features.features_from_address(address, tags={'building':True}, dist=dist)
buildings = osmnx.features.features_from_point(point, tags={'building':True}, dist=dist)
/Users/mtwente/anaconda3/envs/maxvogt/lib/python3.11/site-packages/shapely/set_operations.py:340: RuntimeWarning: invalid value encountered in union
  return lib.union(a, b, **kwargs)
/Users/mtwente/anaconda3/envs/maxvogt/lib/python3.11/site-packages/shapely/predicates.py:798: RuntimeWarning: invalid value encountered in intersects
  return lib.intersects(a, b, **kwargs)
/Users/mtwente/anaconda3/envs/maxvogt/lib/python3.11/site-packages/shapely/predicates.py:798: RuntimeWarning: invalid value encountered in intersects
  return lib.intersects(a, b, **kwargs)
/Users/mtwente/anaconda3/envs/maxvogt/lib/python3.11/site-packages/shapely/predicates.py:798: RuntimeWarning: invalid value encountered in intersects
  return lib.intersects(a, b, **kwargs)
/Users/mtwente/anaconda3/envs/maxvogt/lib/python3.11/site-packages/shapely/predicates.py:798: RuntimeWarning: invalid value encountered in intersects
  return lib.intersects(a, b, **kwargs)
buildings
geometry addr:housenumber addr:street amenity parking access layer addr:city addr:postcode wheelchair ... unisex wheelchair:description recycling:cartons recycling:clothes recycling:electrical_appliances recycling:paper recycling:scrap_metal recycling:wood recycling_type shelter_type
element_type osmid
way 23646201 POLYGON ((8.56476 47.29619, 8.56473 47.29617, ... 21 Bahnhofstrasse NaN NaN NaN NaN Thalwil 8800 NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
28573765 POLYGON ((8.58450 47.31476, 8.58462 47.31461, ... 33 Untere Heslibachstrasse townhall NaN NaN NaN Küsnacht 8700 NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
31316747 POLYGON ((8.56483 47.29642, 8.56498 47.29628, ... 20 Bahnhofstrasse NaN NaN NaN NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
31755378 POLYGON ((8.56547 47.29537, 8.56584 47.29504, ... 15 Bahnhofstrasse NaN NaN NaN NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
31755382 POLYGON ((8.56575 47.29553, 8.56589 47.29538, ... 16 Bahnhofstrasse NaN NaN NaN NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
1299373872 POLYGON ((8.59891 47.29748, 8.59895 47.29742, ... NaN NaN shelter NaN NaN NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN public_transport
1299373875 POLYGON ((8.59849 47.29739, 8.59863 47.29743, ... NaN NaN NaN NaN NaN NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
1299373877 POLYGON ((8.59862 47.29761, 8.59868 47.29762, ... NaN NaN shelter NaN NaN NaN NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN public_transport
1299378652 POLYGON ((8.59160 47.30573, 8.59167 47.30563, ... NaN NaN NaN NaN NaN 1 NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
1305367769 POLYGON ((8.58108 47.31894, 8.58118 47.31896, ... NaN NaN NaN NaN NaN 1 NaN NaN NaN ... NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN

4363 rows × 102 columns

buildings.geom_type.value_counts()
Polygon    4363
Name: count, dtype: int64
buildings = buildings[buildings.geom_type == "Polygon"].reset_index(drop=True)
buildings = buildings[["geometry"]].to_crs(local_crs)
buildings["uID"] = range(len(buildings))
buildings
geometry uID
0 POLYGON ((2685173.904 1238977.426, 2685171.899... 0
1 POLYGON ((2686637.058 1241064.069, 2686646.446... 1
2 POLYGON ((2685178.913 1239003.628, 2685190.829... 2
3 POLYGON ((2685229.206 1238887.977, 2685257.681... 3
4 POLYGON ((2685249.952 1238906.099, 2685261.103... 4
... ... ...
4358 POLYGON ((2687754.918 1239159.357, 2687757.411... 4358
4359 POLYGON ((2687723.328 1239148.749, 2687733.252... 4359
4360 POLYGON ((2687732.266 1239172.499, 2687736.995... 4360
4361 POLYGON ((2687187.977 1240067.743, 2687193.459... 4361
4362 POLYGON ((2686371.382 1241524.644, 2686379.300... 4362

4363 rows × 2 columns

Streets

In comparison to Martin Fleischmann’s workshop, I here set truncate_by_edge to be True. Depending on the OSM mapping quality of the place at hand, it might also be necessary to change the value of network_type. See geopandas documentation for options.

#osm_graph = osmnx.graph_from_place(place, network_type='drive', truncate_by_edge=True)
#osm_graph = osmnx.graph_from_address(place, network_type='drive', truncate_by_edge=True, dist=dist)
osm_graph = osmnx.graph_from_point(point, network_type='drive_service', truncate_by_edge=True, dist=dist)
/Users/mtwente/anaconda3/envs/maxvogt/lib/python3.11/site-packages/osmnx/graph.py:191: FutureWarning: The expected order of coordinates in `bbox` will change in the v2.0.0 release to `(left, bottom, right, top)`.
  G = graph_from_bbox(
/Users/mtwente/anaconda3/envs/maxvogt/lib/python3.11/site-packages/shapely/predicates.py:798: RuntimeWarning: invalid value encountered in intersects
  return lib.intersects(a, b, **kwargs)
/Users/mtwente/anaconda3/envs/maxvogt/lib/python3.11/site-packages/shapely/set_operations.py:340: RuntimeWarning: invalid value encountered in union
  return lib.union(a, b, **kwargs)
/Users/mtwente/anaconda3/envs/maxvogt/lib/python3.11/site-packages/shapely/set_operations.py:340: RuntimeWarning: invalid value encountered in union
  return lib.union(a, b, **kwargs)
osm_graph = osmnx.projection.project_graph(osm_graph, to_crs=local_crs)
streets = osmnx.graph_to_gdfs(
    osm_graph, 
    nodes=False, 
    edges=True,
    node_geometry=False, 
    fill_edge_geometry=True
)
streets.explore()
Make this Notebook Trusted to load map: File -> Trust Notebook
streets
osmid ref name highway maxspeed oneway reversed length geometry lanes service access bridge junction tunnel
u v key
2114177 469665306 0 95160707 17 Seestrasse primary 60 False True 29.118 LINESTRING (2687701.620 1239166.553, 2687712.7... NaN NaN NaN NaN NaN NaN
7602344458 0 813850131 NaN Winkelstrasse service NaN False False 21.408 LINESTRING (2687701.620 1239166.553, 2687681.3... NaN NaN NaN NaN NaN NaN
1833883483 0 95160707 17 Seestrasse primary 60 False False 233.526 LINESTRING (2687701.620 1239166.553, 2687685.6... NaN NaN NaN NaN NaN NaN
469665306 2114177 0 95160707 17 Seestrasse primary 60 False False 29.118 LINESTRING (2687712.706 1239139.619, 2687701.6... NaN NaN NaN NaN NaN NaN
469665338 0 [813850130, 39208990] NaN [Seeweg, Winkelstrasse] residential NaN False False 170.071 LINESTRING (2687712.706 1239139.619, 2687673.6... NaN NaN NaN NaN NaN NaN
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
8942675043 8942675042 0 966590343 NaN NaN service NaN False False 12.543 LINESTRING (2686384.821 1240934.065, 2686391.5... NaN driveway NaN NaN NaN NaN
9356277420 9356277418 0 1014126563 NaN NaN service NaN False True 54.146 LINESTRING (2687189.873 1240347.863, 2687187.1... NaN driveway NaN NaN NaN NaN
10204016687 10204016682 0 1115722503 NaN NaN service NaN False True 63.485 LINESTRING (2686857.333 1241295.033, 2686861.8... NaN NaN NaN NaN NaN NaN
10536484031 10536484032 0 1130572951 NaN NaN service NaN False False 20.476 LINESTRING (2686119.998 1241668.618, 2686128.5... NaN driveway NaN NaN NaN NaN
10810325757 10810325752 0 1162375964 NaN NaN service NaN False False 36.690 LINESTRING (2687639.273 1241235.237, 2687647.9... NaN NaN NaN NaN NaN NaN

2160 rows × 15 columns

streets = momepy.remove_false_nodes(streets)
streets = streets[["geometry"]]
streets["nID"] = range(len(streets))
streets
geometry nID
0 LINESTRING (2687701.620 1239166.553, 2687712.7... 0
1 LINESTRING (2687701.620 1239166.553, 2687685.6... 1
2 LINESTRING (2687712.706 1239139.619, 2687701.6... 2
3 LINESTRING (2687712.706 1239139.619, 2687718.4... 3
4 LINESTRING (2687552.024 1239345.242, 2687561.2... 4
... ... ...
1807 LINESTRING (2688839.282 1240250.976, 2688844.6... 1807
1808 LINESTRING (2689287.152 1241880.389, 2689282.7... 1808
1809 LINESTRING (2689124.751 1238811.304, 2689137.4... 1809
1810 LINESTRING (2688909.001 1239589.492, 2688910.4... 1810
1811 LINESTRING (2689671.698 1241372.801, 2689659.9... 1811

1812 rows × 2 columns

Generated data

Tessellation

We can generate a spatail unit using Voronoi tessellation with given building footprints.

limit = momepy.buffered_limit(buildings, 100)

tessellation = momepy.Tessellation(buildings, "uID", limit, verbose=False, segment=1)
tessellation = tessellation.tessellation
/Users/mtwente/anaconda3/envs/maxvogt/lib/python3.11/site-packages/shapely/constructive.py:181: RuntimeWarning: invalid value encountered in buffer
  return lib.buffer(
/Users/mtwente/anaconda3/envs/maxvogt/lib/python3.11/site-packages/shapely/set_operations.py:426: RuntimeWarning: invalid value encountered in unary_union
  return lib.unary_union(collections, **kwargs)
/var/folders/vs/1r1r5k5x12v5dj0w9z79tj3c0000gn/T/ipykernel_46992/1328706492.py:3: FutureWarning: Class based API like `momepy.Tessellation` is deprecated. Replace it with `momepy.morphological_tessellation` or `momepy.enclosed_tessellation` to use functional API instead or pin momepy version <1.0. Class-based API will be removed in 1.0. 
  tessellation = momepy.Tessellation(buildings, "uID", limit, verbose=False, segment=1)
/Users/mtwente/anaconda3/envs/maxvogt/lib/python3.11/site-packages/shapely/constructive.py:181: RuntimeWarning: invalid value encountered in buffer
  return lib.buffer(
/Users/mtwente/anaconda3/envs/maxvogt/lib/python3.11/site-packages/shapely/set_operations.py:133: RuntimeWarning: invalid value encountered in intersection
  return lib.intersection(a, b, **kwargs)
/var/folders/vs/1r1r5k5x12v5dj0w9z79tj3c0000gn/T/ipykernel_46992/1328706492.py:3: UserWarning: Tessellation contains MultiPolygon elements. Initial objects should  be edited. `unique_id` of affected elements: [3337, 3339, 3341, 3340, 1362, 1443].
  tessellation = momepy.Tessellation(buildings, "uID", limit, verbose=False, segment=1)

Measure

Measure individual morphometric characters.

Dimensions

buildings["area"] = buildings.area
tessellation["area"] = tessellation.area
streets["length"] = streets.length

Shape

buildings['eri'] = momepy.EquivalentRectangularIndex(buildings).series
/var/folders/vs/1r1r5k5x12v5dj0w9z79tj3c0000gn/T/ipykernel_46992/1902756608.py:1: FutureWarning: Class based API like `momepy.EquivalentRectangularIndex` is deprecated. Replace it with `momepy.equivalent_rectangular_index` to use functional API instead or pin momepy version <1.0. Class-based API will be removed in 1.0. 
  buildings['eri'] = momepy.EquivalentRectangularIndex(buildings).series
/Users/mtwente/anaconda3/envs/maxvogt/lib/python3.11/site-packages/pandas/core/arraylike.py:492: RuntimeWarning: invalid value encountered in oriented_envelope
  return getattr(ufunc, method)(*new_inputs, **kwargs)
buildings['elongation'] = momepy.Elongation(buildings).series
/var/folders/vs/1r1r5k5x12v5dj0w9z79tj3c0000gn/T/ipykernel_46992/3723954296.py:1: FutureWarning: Class based API like `momepy.Elongation` is deprecated. Replace it with `momepy.elongation` to use functional API instead or pin momepy version <1.0. Class-based API will be removed in 1.0. 
  buildings['elongation'] = momepy.Elongation(buildings).series
/Users/mtwente/anaconda3/envs/maxvogt/lib/python3.11/site-packages/pandas/core/arraylike.py:492: RuntimeWarning: invalid value encountered in oriented_envelope
  return getattr(ufunc, method)(*new_inputs, **kwargs)
tessellation['convexity'] = momepy.Convexity(tessellation).series
/var/folders/vs/1r1r5k5x12v5dj0w9z79tj3c0000gn/T/ipykernel_46992/2335586521.py:1: FutureWarning: Class based API like `momepy.Convexity` is deprecated. Replace it with `momepy.convexity` to use functional API instead or pin momepy version <1.0. Class-based API will be removed in 1.0. 
  tessellation['convexity'] = momepy.Convexity(tessellation).series
streets["linearity"] = momepy.Linearity(streets).series
/var/folders/vs/1r1r5k5x12v5dj0w9z79tj3c0000gn/T/ipykernel_46992/3802614628.py:1: FutureWarning: Class based API like `momepy.Linearity` is deprecated. Replace it with `momepy.linearity` to use functional API instead or pin momepy version <1.0. Class-based API will be removed in 1.0. 
  streets["linearity"] = momepy.Linearity(streets).series
fig, ax = plt.subplots(1, 2, figsize=(24, 12))

buildings.plot("eri", ax=ax[0], scheme="natural_breaks", legend=True)
buildings.plot("elongation", ax=ax[1], scheme="natural_breaks", legend=True)

ax[0].set_title('Building Equivalent\nRectangular Index', fontsize=20)
ax[1].set_title('Building Elongation', fontsize=20)

ax[0].set_axis_off()
ax[1].set_axis_off()

plt.savefig('../results/Erlenbach/eri_and_elongation.svg')  
/Users/mtwente/anaconda3/envs/maxvogt/lib/python3.11/site-packages/IPython/core/pylabtools.py:77: DeprecationWarning: backend2gui is deprecated since IPython 8.24, backends are managed in matplotlib and can be externally registered.
  warnings.warn(
/Users/mtwente/anaconda3/envs/maxvogt/lib/python3.11/site-packages/IPython/core/pylabtools.py:77: DeprecationWarning: backend2gui is deprecated since IPython 8.24, backends are managed in matplotlib and can be externally registered.
  warnings.warn(

fig, ax = plt.subplots(1, 2, figsize=(24, 12))

tessellation.plot("convexity", ax=ax[0], scheme="natural_breaks", legend=True)
streets.plot("linearity", ax=ax[1], scheme="natural_breaks", legend=True)

ax[0].set_title('Convexity', fontsize=20)
ax[1].set_title('Linearity', fontsize=20)

ax[0].set_axis_off()
ax[1].set_axis_off()

plt.savefig('../results/Erlenbach/convexity_and_linearity.svg')

Spatial distribution

buildings["shared_walls"] = momepy.SharedWallsRatio(buildings).series
/Users/mtwente/anaconda3/envs/maxvogt/lib/python3.11/site-packages/momepy/distribution.py:214: FutureWarning: Class based API like `momepy.SharedWalls` or `momepy.SharedWallsRatio` is deprecated. Replace it with `momepy.shared_walls` or explicitly computing `momepy.shared_walls / gdf.length` respectively to use functional API instead or pin momepy version <1.0. Class-based API will be removed in 1.0. 
  super().__init__(gdf)
/Users/mtwente/anaconda3/envs/maxvogt/lib/python3.11/site-packages/shapely/set_operations.py:133: RuntimeWarning: invalid value encountered in intersection
  return lib.intersection(a, b, **kwargs)
buildings.plot("shared_walls", figsize=(12, 12), scheme="natural_breaks", legend=True).set_axis_off()

plt.savefig('../results/Erlenbach/sharedwalls.svg')

Generate spatial weights matrix using libpysal.

queen_1 = libpysal.weights.contiguity.Queen.from_dataframe(tessellation, ids="uID", silence_warnings=True)
tessellation["neighbors"] = momepy.Neighbors(tessellation, queen_1, "uID", weighted=True, verbose=False).series
tessellation["covered_area"] = momepy.CoveredArea(tessellation, queen_1, "uID", verbose=False).series

with warnings.catch_warnings():
    warnings.simplefilter("ignore")
    
    buildings["neighbor_distance"] = momepy.NeighborDistance(buildings, queen_1, "uID", verbose=False).series
/var/folders/vs/1r1r5k5x12v5dj0w9z79tj3c0000gn/T/ipykernel_46992/1669107924.py:1: FutureWarning: Class based API like `momepy.Neighbors` is deprecated. Replace it with `momepy.neighbors` to use functional API instead or pin momepy version <1.0. Class-based API will be removed in 1.0. 
  tessellation["neighbors"] = momepy.Neighbors(tessellation, queen_1, "uID", weighted=True, verbose=False).series
/var/folders/vs/1r1r5k5x12v5dj0w9z79tj3c0000gn/T/ipykernel_46992/1669107924.py:2: FutureWarning: `momepy.CoveredArea` is deprecated. Replace it with `.describe()` method of libpysal.graph.Graph or pin momepy version <1.0. This class will be removed in 1.0. 
  tessellation["covered_area"] = momepy.CoveredArea(tessellation, queen_1, "uID", verbose=False).series
fig, ax = plt.subplots(1, 2, figsize=(24, 12))

buildings.plot("neighbor_distance", ax=ax[0], scheme="natural_breaks", legend=True)
tessellation.plot("covered_area", ax=ax[1], scheme="natural_breaks", legend=True)

ax[0].set_title('Neighbour Distance', fontsize=20)
ax[1].set_title('Covered Area', fontsize=20)


ax[0].set_axis_off()
ax[1].set_axis_off()

plt.savefig('../results/Erlenbach/neighbourdist_and_coveredarea.svg')  

queen_3 = momepy.sw_high(k=3, weights=queen_1)
buildings_q1 = libpysal.weights.contiguity.Queen.from_dataframe(buildings, silence_warnings=True)

buildings['interbuilding_distance'] = momepy.MeanInterbuildingDistance(buildings, queen_1, 'uID', queen_3, verbose=False).series
buildings['adjacency'] = momepy.BuildingAdjacency(buildings, queen_3, 'uID', buildings_q1, verbose=False).series
/var/folders/vs/1r1r5k5x12v5dj0w9z79tj3c0000gn/T/ipykernel_46992/3221532125.py:1: FutureWarning: `momepy.sw_high` is deprecated. Replace it with .higher_order() method of libpysal.graph.Graph or pin momepy version <1.0. This class will be removed in 1.0. 
  queen_3 = momepy.sw_high(k=3, weights=queen_1)
/var/folders/vs/1r1r5k5x12v5dj0w9z79tj3c0000gn/T/ipykernel_46992/3221532125.py:2: FutureWarning: `use_index` defaults to False but will default to True in future. Set True/False directly to control this behavior and silence this warning
  buildings_q1 = libpysal.weights.contiguity.Queen.from_dataframe(buildings, silence_warnings=True)
/var/folders/vs/1r1r5k5x12v5dj0w9z79tj3c0000gn/T/ipykernel_46992/3221532125.py:4: FutureWarning: Class based API like `momepy.MeanInterbuildingDistance` is deprecated. Replace it with `momepy.mean_interbuilding_distance` to use functional API instead or pin momepy version <1.0. Class-based API will be removed in 1.0. 
  buildings['interbuilding_distance'] = momepy.MeanInterbuildingDistance(buildings, queen_1, 'uID', queen_3, verbose=False).series
/Users/mtwente/anaconda3/envs/maxvogt/lib/python3.11/site-packages/shapely/measurement.py:74: RuntimeWarning: invalid value encountered in distance
  return lib.distance(a, b, **kwargs)
/var/folders/vs/1r1r5k5x12v5dj0w9z79tj3c0000gn/T/ipykernel_46992/3221532125.py:5: FutureWarning: Class based API like `momepy.BuildingAdjacency` is deprecated. Replace it with `momepy.building_adjacency` to use functional API instead or pin momepy version <1.0. Class-based API will be removed in 1.0. 
  buildings['adjacency'] = momepy.BuildingAdjacency(buildings, queen_3, 'uID', buildings_q1, verbose=False).series
#fig, ax = plt.subplots(1, 2, figsize=(24, 12))
#
#buildings.plot("interbuilding_distance", ax=ax[0], scheme="natural_breaks", legend=True)
#buildings.plot("adjacency", ax=ax[1], scheme="natural_breaks", legend=True)
#
#ax[0].set_axis_off()
#ax[1].set_axis_off()
profile = momepy.StreetProfile(streets, buildings)
streets["width"] = profile.w
streets["width_deviation"] = profile.wd
streets["openness"] = profile.o
/var/folders/vs/1r1r5k5x12v5dj0w9z79tj3c0000gn/T/ipykernel_46992/2648127835.py:1: FutureWarning: Class based API like `momepy.StreetProfile` is deprecated. Replace it with `momepy.street_profile` to use functional API instead or pin momepy version <1.0. Class-based API will be removed in 1.0. 
  profile = momepy.StreetProfile(streets, buildings)
/Users/mtwente/anaconda3/envs/maxvogt/lib/python3.11/site-packages/pandas/core/arraylike.py:492: RuntimeWarning: invalid value encountered in intersection
  return getattr(ufunc, method)(*new_inputs, **kwargs)
fig, ax = plt.subplots(1, 3, figsize=(24, 12))

streets.plot("width", ax=ax[0], scheme="natural_breaks", legend=True)
streets.plot("width_deviation", ax=ax[1], scheme="natural_breaks", legend=True)
streets.plot("openness", ax=ax[2], scheme="natural_breaks", legend=True)

ax[0].set_title('Width', fontsize=20)
ax[1].set_title('Width Deviation', fontsize=20)
ax[2].set_title('Openness', fontsize=20)

ax[0].set_axis_off()
ax[1].set_axis_off()
ax[2].set_axis_off()

plt.savefig('../results/Erlenbach/road_network.svg')

Intensity

tessellation['car'] = momepy.AreaRatio(tessellation, buildings, 'area', 'area', 'uID').series
/var/folders/vs/1r1r5k5x12v5dj0w9z79tj3c0000gn/T/ipykernel_46992/1819725789.py:1: FutureWarning: `momepy.AreaRatio` is deprecated. Replace it with a direct division of areas or momepy.describe_agg() or pin momepy version <1.0. This class will be removed in 1.0. 
  tessellation['car'] = momepy.AreaRatio(tessellation, buildings, 'area', 'area', 'uID').series
tessellation.plot("car", figsize=(12, 12), vmin=0, vmax=1, legend=True).set_axis_off()
plt.title("Building/Tessellation Area Ratio")

plt.savefig('../results/Erlenbach/tessellation_ratio.svg')

Connectivity

graph = momepy.gdf_to_nx(streets)
graph = momepy.node_degree(graph)
graph = momepy.closeness_centrality(graph, radius=400, distance="mm_len")
graph = momepy.meshedness(graph, radius=400, distance="mm_len")
nodes, streets = momepy.nx_to_gdf(graph)
fig, ax = plt.subplots(1, 3, figsize=(24, 12))

nodes.plot("degree", ax=ax[0], scheme="natural_breaks", legend=True, markersize=1)
nodes.plot("closeness", ax=ax[1], scheme="natural_breaks", legend=True, markersize=1, legend_kwds={"fmt": "{:.6f}"})
nodes.plot("meshedness", ax=ax[2], scheme="natural_breaks", legend=True, markersize=1)

ax[0].set_axis_off()
ax[1].set_axis_off()
ax[2].set_axis_off()

plt.savefig('../results/Erlenbach/connectivity.svg')  

buildings["nodeID"] = momepy.get_node_id(buildings, nodes, streets, "nodeID", "nID")
/var/folders/vs/1r1r5k5x12v5dj0w9z79tj3c0000gn/T/ipykernel_46992/403706460.py:1: FutureWarning: Class based API like `momepy.get_node_id` is deprecated. Replace it with `momepy.get_nearest_node` to use functional API instead or pin momepy version <1.0. Class-based API will be removed in 1.0. 
  buildings["nodeID"] = momepy.get_node_id(buildings, nodes, streets, "nodeID", "nID")

Link all data together (to tessellation cells or buildings).

tessellation
uID geometry nID area convexity neighbors covered_area car
0 3134 POLYGON ((2688829.623 1237987.800, 2688832.637... 1791.0 19789.953804 0.943741 0.009088 37331.573639 0.078231
1 3235 POLYGON ((2689203.090 1238093.317, 2689202.652... 1802.0 15539.551230 0.992607 0.010454 34446.760704 0.008696
2 3266 POLYGON ((2689153.845 1238089.253, 2689153.425... 1802.0 7278.007503 0.917167 0.010546 29643.654087 0.063659
3 3237 POLYGON ((2689112.133 1238088.000, 2689111.933... 1798.0 4647.593610 0.966758 0.012934 17485.918532 0.031800
4 3274 POLYGON ((2689041.855 1238075.424, 2689041.818... 1798.0 4505.267563 0.987124 0.013113 16500.334713 0.028477
... ... ... ... ... ... ... ... ...
4358 1314 POLYGON ((2687003.304 1241983.900, 2687002.273... 1597.0 770.928778 0.945052 0.042760 8532.906174 0.191557
4359 1346 POLYGON ((2686980.837 1241988.704, 2686980.050... 1597.0 1229.455390 0.873244 0.040315 20119.461617 0.143980
4360 1319 POLYGON ((2686950.978 1242029.662, 2686950.870... 1585.0 3264.755115 0.953183 0.014082 26006.578160 0.066619
4361 1328 POLYGON ((2686997.021 1242040.707, 2686993.745... 1597.0 8574.958447 0.942077 0.012853 27950.589254 0.034091
4362 3738 POLYGON ((2687099.056 1242060.927, 2687099.216... 1611.0 15577.245146 0.926410 0.007921 27637.484607 0.036476

4363 rows × 8 columns

merged = tessellation.merge(buildings.drop(columns=['nID', 'geometry']), on='uID')
merged = merged.merge(streets.drop(columns='geometry'), on='nID', how='left')
merged = merged.merge(nodes.drop(columns='geometry'), on='nodeID', how='left')
merged.columns
Index(['uID', 'geometry', 'nID', 'area_x', 'convexity', 'neighbors',
       'covered_area', 'car', 'area_y', 'eri', 'elongation', 'shared_walls',
       'neighbor_distance', 'interbuilding_distance', 'adjacency', 'nodeID',
       'length', 'linearity', 'width', 'width_deviation', 'openness', 'mm_len',
       'node_start', 'node_end', 'x', 'y', 'degree', 'closeness',
       'meshedness'],
      dtype='object')

Understanding the context

Measure first, second and third quartile of distribution of values within an area around each building.

percentiles = []
for column in merged.columns.drop(["uID", "nodeID", "nID", 'mm_len', 'node_start', 'node_end', "geometry"]):
    perc = momepy.Percentiles(merged, column, queen_3, "uID", verbose=False).frame
    perc.columns = [f"{column}_" + str(x) for x in perc.columns]
    percentiles.append(perc)
/var/folders/vs/1r1r5k5x12v5dj0w9z79tj3c0000gn/T/ipykernel_46992/962013058.py:3: FutureWarning: Class based API like `momepy.Percentiles` is deprecated. Replace it with `momepy.percentile` to use functional API instead or pin momepy version <1.0. Class-based API will be removed in 1.0. 
  perc = momepy.Percentiles(merged, column, queen_3, "uID", verbose=False).frame
/var/folders/vs/1r1r5k5x12v5dj0w9z79tj3c0000gn/T/ipykernel_46992/962013058.py:3: FutureWarning: Class based API like `momepy.Percentiles` is deprecated. Replace it with `momepy.percentile` to use functional API instead or pin momepy version <1.0. Class-based API will be removed in 1.0. 
  perc = momepy.Percentiles(merged, column, queen_3, "uID", verbose=False).frame
/var/folders/vs/1r1r5k5x12v5dj0w9z79tj3c0000gn/T/ipykernel_46992/962013058.py:3: FutureWarning: Class based API like `momepy.Percentiles` is deprecated. Replace it with `momepy.percentile` to use functional API instead or pin momepy version <1.0. Class-based API will be removed in 1.0. 
  perc = momepy.Percentiles(merged, column, queen_3, "uID", verbose=False).frame
/var/folders/vs/1r1r5k5x12v5dj0w9z79tj3c0000gn/T/ipykernel_46992/962013058.py:3: FutureWarning: Class based API like `momepy.Percentiles` is deprecated. Replace it with `momepy.percentile` to use functional API instead or pin momepy version <1.0. Class-based API will be removed in 1.0. 
  perc = momepy.Percentiles(merged, column, queen_3, "uID", verbose=False).frame
/var/folders/vs/1r1r5k5x12v5dj0w9z79tj3c0000gn/T/ipykernel_46992/962013058.py:3: FutureWarning: Class based API like `momepy.Percentiles` is deprecated. Replace it with `momepy.percentile` to use functional API instead or pin momepy version <1.0. Class-based API will be removed in 1.0. 
  perc = momepy.Percentiles(merged, column, queen_3, "uID", verbose=False).frame
/var/folders/vs/1r1r5k5x12v5dj0w9z79tj3c0000gn/T/ipykernel_46992/962013058.py:3: FutureWarning: Class based API like `momepy.Percentiles` is deprecated. Replace it with `momepy.percentile` to use functional API instead or pin momepy version <1.0. Class-based API will be removed in 1.0. 
  perc = momepy.Percentiles(merged, column, queen_3, "uID", verbose=False).frame
/var/folders/vs/1r1r5k5x12v5dj0w9z79tj3c0000gn/T/ipykernel_46992/962013058.py:3: FutureWarning: Class based API like `momepy.Percentiles` is deprecated. Replace it with `momepy.percentile` to use functional API instead or pin momepy version <1.0. Class-based API will be removed in 1.0. 
  perc = momepy.Percentiles(merged, column, queen_3, "uID", verbose=False).frame
/var/folders/vs/1r1r5k5x12v5dj0w9z79tj3c0000gn/T/ipykernel_46992/962013058.py:3: FutureWarning: Class based API like `momepy.Percentiles` is deprecated. Replace it with `momepy.percentile` to use functional API instead or pin momepy version <1.0. Class-based API will be removed in 1.0. 
  perc = momepy.Percentiles(merged, column, queen_3, "uID", verbose=False).frame
/var/folders/vs/1r1r5k5x12v5dj0w9z79tj3c0000gn/T/ipykernel_46992/962013058.py:3: FutureWarning: Class based API like `momepy.Percentiles` is deprecated. Replace it with `momepy.percentile` to use functional API instead or pin momepy version <1.0. Class-based API will be removed in 1.0. 
  perc = momepy.Percentiles(merged, column, queen_3, "uID", verbose=False).frame
/var/folders/vs/1r1r5k5x12v5dj0w9z79tj3c0000gn/T/ipykernel_46992/962013058.py:3: FutureWarning: Class based API like `momepy.Percentiles` is deprecated. Replace it with `momepy.percentile` to use functional API instead or pin momepy version <1.0. Class-based API will be removed in 1.0. 
  perc = momepy.Percentiles(merged, column, queen_3, "uID", verbose=False).frame
/var/folders/vs/1r1r5k5x12v5dj0w9z79tj3c0000gn/T/ipykernel_46992/962013058.py:3: FutureWarning: Class based API like `momepy.Percentiles` is deprecated. Replace it with `momepy.percentile` to use functional API instead or pin momepy version <1.0. Class-based API will be removed in 1.0. 
  perc = momepy.Percentiles(merged, column, queen_3, "uID", verbose=False).frame
/var/folders/vs/1r1r5k5x12v5dj0w9z79tj3c0000gn/T/ipykernel_46992/962013058.py:3: FutureWarning: Class based API like `momepy.Percentiles` is deprecated. Replace it with `momepy.percentile` to use functional API instead or pin momepy version <1.0. Class-based API will be removed in 1.0. 
  perc = momepy.Percentiles(merged, column, queen_3, "uID", verbose=False).frame
/var/folders/vs/1r1r5k5x12v5dj0w9z79tj3c0000gn/T/ipykernel_46992/962013058.py:3: FutureWarning: Class based API like `momepy.Percentiles` is deprecated. Replace it with `momepy.percentile` to use functional API instead or pin momepy version <1.0. Class-based API will be removed in 1.0. 
  perc = momepy.Percentiles(merged, column, queen_3, "uID", verbose=False).frame
/var/folders/vs/1r1r5k5x12v5dj0w9z79tj3c0000gn/T/ipykernel_46992/962013058.py:3: FutureWarning: Class based API like `momepy.Percentiles` is deprecated. Replace it with `momepy.percentile` to use functional API instead or pin momepy version <1.0. Class-based API will be removed in 1.0. 
  perc = momepy.Percentiles(merged, column, queen_3, "uID", verbose=False).frame
/var/folders/vs/1r1r5k5x12v5dj0w9z79tj3c0000gn/T/ipykernel_46992/962013058.py:3: FutureWarning: Class based API like `momepy.Percentiles` is deprecated. Replace it with `momepy.percentile` to use functional API instead or pin momepy version <1.0. Class-based API will be removed in 1.0. 
  perc = momepy.Percentiles(merged, column, queen_3, "uID", verbose=False).frame
/var/folders/vs/1r1r5k5x12v5dj0w9z79tj3c0000gn/T/ipykernel_46992/962013058.py:3: FutureWarning: Class based API like `momepy.Percentiles` is deprecated. Replace it with `momepy.percentile` to use functional API instead or pin momepy version <1.0. Class-based API will be removed in 1.0. 
  perc = momepy.Percentiles(merged, column, queen_3, "uID", verbose=False).frame
/var/folders/vs/1r1r5k5x12v5dj0w9z79tj3c0000gn/T/ipykernel_46992/962013058.py:3: FutureWarning: Class based API like `momepy.Percentiles` is deprecated. Replace it with `momepy.percentile` to use functional API instead or pin momepy version <1.0. Class-based API will be removed in 1.0. 
  perc = momepy.Percentiles(merged, column, queen_3, "uID", verbose=False).frame
/var/folders/vs/1r1r5k5x12v5dj0w9z79tj3c0000gn/T/ipykernel_46992/962013058.py:3: FutureWarning: Class based API like `momepy.Percentiles` is deprecated. Replace it with `momepy.percentile` to use functional API instead or pin momepy version <1.0. Class-based API will be removed in 1.0. 
  perc = momepy.Percentiles(merged, column, queen_3, "uID", verbose=False).frame
/var/folders/vs/1r1r5k5x12v5dj0w9z79tj3c0000gn/T/ipykernel_46992/962013058.py:3: FutureWarning: Class based API like `momepy.Percentiles` is deprecated. Replace it with `momepy.percentile` to use functional API instead or pin momepy version <1.0. Class-based API will be removed in 1.0. 
  perc = momepy.Percentiles(merged, column, queen_3, "uID", verbose=False).frame
/var/folders/vs/1r1r5k5x12v5dj0w9z79tj3c0000gn/T/ipykernel_46992/962013058.py:3: FutureWarning: Class based API like `momepy.Percentiles` is deprecated. Replace it with `momepy.percentile` to use functional API instead or pin momepy version <1.0. Class-based API will be removed in 1.0. 
  perc = momepy.Percentiles(merged, column, queen_3, "uID", verbose=False).frame
/var/folders/vs/1r1r5k5x12v5dj0w9z79tj3c0000gn/T/ipykernel_46992/962013058.py:3: FutureWarning: Class based API like `momepy.Percentiles` is deprecated. Replace it with `momepy.percentile` to use functional API instead or pin momepy version <1.0. Class-based API will be removed in 1.0. 
  perc = momepy.Percentiles(merged, column, queen_3, "uID", verbose=False).frame
/var/folders/vs/1r1r5k5x12v5dj0w9z79tj3c0000gn/T/ipykernel_46992/962013058.py:3: FutureWarning: Class based API like `momepy.Percentiles` is deprecated. Replace it with `momepy.percentile` to use functional API instead or pin momepy version <1.0. Class-based API will be removed in 1.0. 
  perc = momepy.Percentiles(merged, column, queen_3, "uID", verbose=False).frame
percentiles_joined = pandas.concat(percentiles, axis=1)
percentiles_joined
area_x_25 area_x_50 area_x_75 convexity_25 convexity_50 convexity_75 neighbors_25 neighbors_50 neighbors_75 covered_area_25 ... y_75 degree_25 degree_50 degree_75 closeness_25 closeness_50 closeness_75 meshedness_25 meshedness_50 meshedness_75
0 1138.012395 1870.719366 4794.388658 0.942479 0.962164 0.974379 0.017494 0.032391 0.039986 12329.610123 ... 1.238251e+06 1.0 1.0 4.0 0.000005 0.000006 0.000029 0.000000 0.462185 1.000000
1 912.287838 1324.257689 3631.333462 0.938405 0.955914 0.972011 0.018805 0.038831 0.049242 12551.404477 ... 1.238274e+06 1.0 1.0 4.0 0.000044 0.000050 0.000064 0.268293 0.304348 0.333333
2 852.974597 1304.932154 4068.300513 0.938417 0.959628 0.970982 0.017407 0.039152 0.047123 10796.059015 ... 1.238214e+06 1.0 1.0 1.0 0.000044 0.000044 0.000052 0.296296 0.333333 0.333333
3 844.796549 1138.105912 4068.300513 0.938417 0.962806 0.973349 0.018216 0.040810 0.053376 7852.309954 ... 1.238168e+06 1.0 1.0 1.0 0.000005 0.000044 0.000044 0.000000 0.333333 0.333333
4 866.633922 1158.043873 3044.249730 0.949214 0.962751 0.972805 0.021715 0.040583 0.048357 9177.682742 ... 1.238037e+06 1.0 1.0 1.0 0.000005 0.000006 0.000006 0.000000 0.166667 0.666667
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
4358 1258.199607 1661.545491 2266.896633 0.930152 0.946325 0.965933 0.029686 0.036613 0.042401 9865.651107 ... 1.242028e+06 1.0 5.0 5.5 0.000075 0.000123 0.000136 0.418182 0.428571 0.431373
4359 1258.199607 1853.848992 2343.252672 0.937895 0.951443 0.967834 0.027272 0.035514 0.042401 10217.356669 ... 1.242028e+06 1.0 4.5 5.0 0.000075 0.000102 0.000134 0.418182 0.428571 0.432667
4360 1279.901646 1976.668136 2781.339587 0.933947 0.945421 0.966140 0.013732 0.035277 0.042401 11056.133159 ... 1.242028e+06 1.0 1.0 5.0 0.000056 0.000075 0.000123 0.418182 0.428571 0.448276
4361 1308.645863 1941.034346 2781.339587 0.930152 0.951443 0.966140 0.013732 0.034213 0.041067 10173.174499 ... 1.242028e+06 1.0 1.0 5.0 0.000065 0.000088 0.000127 0.418182 0.428571 0.445068
4362 1127.817342 1424.690469 2244.426001 0.894416 0.936052 0.958186 0.021509 0.035550 0.041067 8922.859123 ... 1.242028e+06 1.0 5.0 6.0 0.000077 0.000103 0.000127 0.400000 0.416408 0.428571

4363 rows × 66 columns

fig, ax = plt.subplots(1, 2, figsize=(24, 12))

tessellation.plot("convexity", ax=ax[0], scheme="natural_breaks", legend=True)
merged.plot(percentiles_joined['convexity_50'].values, ax=ax[1], scheme="natural_breaks", legend=True)

ax[0].set_title('Convexity', fontsize=20)
ax[1].set_title('Convexity_50', fontsize=20)

ax[0].set_axis_off()
ax[1].set_axis_off()

plt.savefig('../results/Erlenbach/convexity.svg')

Clustering

Standardize values before clustering.

standardized = (percentiles_joined - percentiles_joined.mean()) / percentiles_joined.std()
standardized
area_x_25 area_x_50 area_x_75 convexity_25 convexity_50 convexity_75 neighbors_25 neighbors_50 neighbors_75 covered_area_25 ... y_75 degree_25 degree_50 degree_75 closeness_25 closeness_50 closeness_75 meshedness_25 meshedness_50 meshedness_75
0 0.157751 0.304221 1.187585 0.862473 0.870777 0.303628 -1.825162 -1.155740 -1.059137 0.874604 ... -1.836746 -0.977 -2.968876 -1.968592 -1.471781 -1.639327 -1.475981 -3.966363 -0.078288 4.109512
1 -0.025760 -0.052042 0.640771 0.609642 0.190264 -0.093325 -1.675229 -0.361398 -0.052160 0.922370 ... -1.816922 -0.977 -2.968876 -1.968592 -0.818382 -1.011104 -1.027846 -1.488803 -1.557979 -1.524357
2 -0.073981 -0.064641 0.846213 0.610381 0.594625 -0.265926 -1.835181 -0.321829 -0.282702 0.544334 ... -1.869982 -0.977 -2.968876 -5.392496 -0.818382 -1.086965 -1.182698 -1.230203 -1.286246 -1.524357
3 -0.080629 -0.173402 0.846213 0.610381 0.940615 0.130870 -1.742562 -0.117307 0.397582 -0.089641 ... -1.910532 -0.977 -2.968876 -5.392496 -1.471781 -1.086965 -1.281266 -3.966363 -1.286246 -1.524357
4 -0.062876 -0.160404 0.364752 1.280374 0.934586 0.039748 -1.342325 -0.145327 -0.148436 0.195795 ... -2.027066 -0.977 -2.968876 -5.392496 -1.471781 -1.643678 -1.773867 -3.966363 -2.848713 1.292578
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
4358 0.255461 0.167851 -0.000723 0.097490 -0.853641 -1.112512 -0.430429 -0.634927 -0.796447 0.343958 ... 1.518616 -0.977 0.107566 -0.256640 -0.310440 0.046045 -0.114111 -0.104646 -0.393407 -0.695846
4359 0.255461 0.293223 0.035176 0.577998 -0.296476 -0.793689 -0.706585 -0.770514 -0.796447 0.419703 ... 1.518616 -0.977 -0.276989 -0.827291 -0.310440 -0.264948 -0.137974 -0.104646 -0.393407 -0.684904
4360 0.273104 0.373294 0.241144 0.333010 -0.952107 -1.077856 -2.255579 -0.799804 -0.796447 0.600344 ... 1.518616 -0.977 -2.968876 -0.827291 -0.619077 -0.650804 -0.270839 -0.104646 -0.393407 -0.553000
4361 0.296473 0.350063 0.241144 0.097490 -0.296476 -1.077856 -2.255579 -0.930974 -0.941474 0.410187 ... 1.518616 -0.977 -2.968876 -0.827291 -0.466483 -0.457402 -0.224959 -0.104646 -0.393407 -0.580107
4362 0.149462 0.013435 -0.011288 -2.120159 -1.972012 -2.411490 -1.365830 -0.766097 -0.941474 0.140915 ... 1.518616 -0.977 0.107566 0.314010 -0.278265 -0.241244 -0.224959 -0.272547 -0.507437 -0.719518

4363 rows × 66 columns

How many clusters?

cgram = Clustergram(range(1, 12), n_init=10, random_state=0)
cgram.fit(standardized.fillna(0))
K=1 skipped. Mean computed from data directly.
K=2 fitted in 0.071 seconds.
K=3 fitted in 0.053 seconds.
K=4 fitted in 0.065 seconds.
K=5 fitted in 0.149 seconds.
K=6 fitted in 0.092 seconds.
K=7 fitted in 0.110 seconds.
K=8 fitted in 0.099 seconds.
K=9 fitted in 0.126 seconds.
K=10 fitted in 0.111 seconds.
K=11 fitted in 0.119 seconds.
Clustergram(k_range=range(1, 12), backend='sklearn', method='kmeans', kwargs={'n_init': 10, 'random_state': 0})
show(cgram.bokeh())
cgram.labels.head()
1 2 3 4 5 6 7 8 9 10 11
0 0 0 2 0 3 1 0 1 0 7 10
1 0 0 2 0 3 1 0 6 0 7 10
2 0 0 2 0 3 1 0 6 0 7 10
3 0 0 2 0 3 1 0 6 0 7 10
4 0 0 2 0 3 1 0 1 0 7 10 
merged["cluster"] = cgram.labels[8].values
urban_types = buildings[["geometry", "uID"]].merge(merged[["uID", "cluster"]], on="uID")
urban_types.explore("cluster", categorical=True, prefer_canvas=True, tiles="CartoDB Positron", tooltip=False)
Make this Notebook Trusted to load map: File -> Trust Notebook

Export

Save cluster output geodata to geojson file:

urban_types.to_file('../data/raw/Erlenbach.geojson', driver='GeoJSON')

Utilities

Once you run the cell below, it hides all cells below the currently active one (aka presentation mode).

%%html
<style>
.jp-Cell.jp-mod-selected ~ .jp-Cell {
    display: none;
}
</style>

Run this cell to turn the presentation mode off.

%%html
<style>
.jp-Cell.jp-mod-selected ~ .jp-Cell {
    display: block;
}
</style>
Back to top